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Many ground state studies of 4 He using a shadow wave function with an in- 
verse fifth power McMillan particle-particle correlation function have yielded 

•4— > ' 

radial distribution functions with misplaced peaks. It has been conjectured 

S' 

that this is due to the specific choice of the McMillan correlation function. 
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However, beyond the use of fully optimized two-particle correlation functions, 
there has been little study of simple alternatives that can correct this de- 



y\ ' feet. In this work we show that the remedy is surprisingly simple. When a 



shadow wavefunction with an inverse seventh power particle-particle correla- 
tion function is used to study 4 He, it gives a correctly peaked radial distribu- 
tion function, lowers the energy at all liquid and solid densities, and produces 

an excellent equation of state. 
05.30.Fk, 21.65.+f, 67.55.-s 



Typeset using REVTgX 



The ground state properties of liquid and solid 4 He have been studied extensively over 
the years by a variety of many-body techniques, ranging from hypernetted-chain (HNC) 
theories'™ and variational Monte Carlo (VMC) methods^ to "exact" Monte Carlo meth- 
ods (EMC) such as Green's Function Monte Carlcii (GFMC) and Diffusion Monte Carlo@ 
(DMC). The advance of EMC have seemingly oblivated the need for the "brutish" and 
bias-laiden method of VMC. However, the introduction of the shadow wavefunction by Vi- 
tiello et aM has added new subtlety and refinement to this approach. It has been shown 
that shadow wavefunctions can describe both the liquid and solid phase of 4 He with ex- 
cellent energy, while simultaneously maintaining translational and Bose sysmmetry. Since 
any VMC calculation is an order of magnitude less computationally demanding than cor- 
responding EMC calculations, the method of shadow wavefunctions remained economically 
and conceptionally appealing. 

However, it has been noted for some time that shadow wavefunctions with the McMillan 
inverse fifth power particle-particle correlation function do not give a correct radial distri- 
bution function for bulk liquid 4 He. All the peaks are misplaced as in the original McMillan 
calculation. While there have been continued improvements on the form of two-particle cor- 
relation, leading to fully optimized correlations expressible in terms of a basis state, there 
has been no explorations of simpler alternatives to cure this defect. In this work, we shown 
that a simple replacement of the inverse fifth power by that of an inverse seventh power 
significantly improves simultaneously, the ground state energy, the equation of state and the 
radial distribution function. 

The first VMC calculation of the groundstate properties of liquid and solid 4 He was 
carried out by McMillanH who employed a trial wave function with an inverse fifth power 
of the particle separation as the two-body correlation function. In this early study the 
potential between 4 He atoms is taken to be the two-body Lennard- Jones (LJ) potential, 
v(r) = 4e[(a/r) 12 — (o~/r) 6 ], with the DeBoer-Michels parameters, a = 2.556 A and e = 
10.22K . Since then, the two-body HFDHE2 potential of Aziz et affl has superceded the 
LJ potential as the potential of choice for the Helium studies. More recent minor revisions 



of this potentially have added, but have not greatly altered the quality of description of the 
interaction between Helium atoms at low pressure. In this work, we will continue to use the 
HFDHE2 potential to facilitate comparison with existing calculations in the literature. 
The original shadow wavefunction of Vitiello et aM can be written as: 

*(R) = $(R) J dS 6(R, S) $ S (S) , (1) 



with a Gaussian particle-shadow correlation function 

N 



9(R,S) =exp 
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where R = {ri, r 2 , . . . , r^} and S = {si, s 2 , . . . , s^v} represent the set of particle and shadow 
coordinates. The particle-particle correlation function $(R) and the shadow-shadow corre- 
lation function $s(S) are both of the pair-wise product form, $(R) = exp —Y.i>j u ( r ij] 
and $s(S) = exp —J2i>j u s(sij) , where r^ = jr, — Tj\ and s^- = |sj — Sj\. Both $(R) 
and $s(S) are taken as McMillan form of inverse fifth powers (m=5) with u(r) = \{^) m 
and us(s) = (— ) m , where b and b' are two variational parameters. The one-to-one coupling 
constant C between particles and shadows is treated as a variational parameter. We refer 
to this wavefunction as M+MS. 

Optimizing the shadow-shadow pseudopotentiafcU in the form of us{s) = (— ) n with two 
variational parameters b' and n, produces no significant improvements. A better choice 
was found following a suggestion by Reatto et alc3 : Us{s) = r v(as), where v(r) is the 
Aziz HFDHE2 potential and r and a are variational parameters. This scaled Aziz shadow 
correlation, introduced by MacFarland et alcl, which uses a McMillan inverse fifth power 
particle-particle psedopotential (m=5), will be denoted as M+AS. The M+AS shadow wave- 
function substantially lowered the variational energies and improved the description of liquid 
and solid 4 He at all densities. For example, at the GFMC equilibrium density pa 3 = 0.365 
and at freezing density of pa 3 = 0.438, the M+AS energy is about 0.5 K° lower than the 
M+MS energy. 

Both wavef unctions, however, produce a radial distribution function at equilibrium den- 
sity whose main peak is shifted outward by about 0.1 A as compared with the experimental 



value. The same misplacement is also observed at the GFMC freezing density. Such a mis- 
placement can be corrected by optimizing the two-body correlations through the method of 
basis state expansion^. (The peak height is still underestimated, however.) 

In this work, we show that this crucial defect can be simply corrected by a better choice 
of the inverse power (from 5 to 7) in McMillan's form of the particle-particle correlation 
function. 

Without reoptimizing the M+AS wavef unction's shadow parameters, but only varying 
the variational parameter b with m = 7, we obtained lower energies than those of M+AS at 
all liquid and solid densities, in an amount ranging from 0.1 to 0.3K°. This choice of the 
wavefunction in our work, referred to as M7+AS, allows us to improve the quality of the 
shadow wavefunction while retaining the same level of simplicity as before. 

For a system of N Helium atoms interacting via two-body forces only, the Hamiltonian 
has the form 

t2 JV N 

1=1 l>] 

where v (r^) is the Aziz HFDHE2 potential. The expectation value of the Hamiltonian can 
be expressed as 

JdR|^(R)| 2 l ' 

= f dRdSdS' p(R, S, S') E L (R, S, S') . (5) 

The local energy is written as 

£z(R,S,S)- $(R)e(R5S) , (6) 

and does not depend on $^(S) since H acts only upon the variables describing the system 
of real particles. 

The probability p(R, S, S') is given by 



p(R,S,S0 = 

$(R) 2 8(R,S) gg(S) 9(R,SQ $ S (S') 
JdRdS dS' $(R) 2 0(R,S) $ S (S) e(R,S') $ S (S') 



(7) 



To evaluate the expectation value of the Hamiltonian we use the Metropolis Monte Carlo 
algorithmic to sample the probability density p(R, S, S') from the 9N dimensional config- 
uration space of the particles and two sets of shadow coordinates. In these computations, 
the Metropolis steps are subdivided in two parts. In the first, one attempts to move real 
particle coordinates at random inside cubical boxes of side length A. In the second, analo- 
gous attempts to move shadow coordinates are made inside cubical boxes of side length A'. 
After we attempt to move all the shadow coordinates of set {S}, the same is done for those 
in set {S'}. The parameters A and A' were adjusted so that the acceptance ratio for both 
particle and shadow moves was nearly 50%. 

We compute the ground state variational energy, the radial distribution function g(r), 
and the static structure factor S(k). These quantities are spherical averages and have been 
computed for both the real particles and the shadow coordinates. The radial distribution 
function is defined by 

9(r) = ^-J:(S(\r l -r 3 -r\)), (8) 

where the angular brackets denote an average with respect to |\I/(R)| 2 and p is the particle 
density. The static structure factor is obtained from the average 

S(k) = ^<p- k p k > , (9) 

where pk is given by p k = X/^=i ex P( — *kij). By using this procedure S(k) is computed 
for a discrete set of k values where the smaller wave vector compatible with the periodic 
boundary condition of the system is k = 2tt/L ( L is the side of the simulation box ). 

All simulations presented in this work have been done with N = 108 atoms of 4 He 
in a cubic box with periodic boundary conditions. To enforce periodicity, the two-body 
interaction potential v(r) smoothly goes to zero at a cutoff distance, r c = L/2, equal half 
the side of the simulation box. We actually use a slightly modified two-body interaction 
potential v'(r) = v(r) — Av(r) according to the replacement 



v(r) + v (2r c — r) — 2v(r c ), r < r c 
v\r) = { ~ (10) 

0, r > r c . 

A correction AV = (p/2) J d 3 rg(r)Av (r) was then added to the computed potential energy, 
where the radial distribution function g(r) comes from the simulation and is taken equal to 
1 for r > r c . The shadow-shadow pseudopotential us(s) was modified according to the same 
prescription as v(r), while the particle-particle inverse power McMillan pseudopotential u(r) 
and its first two derivatives were slightly modified near the edge of the simulation box in 
order to go smoothly to zero, by using a third degree polynomial fit to the pseudopotential 
near the edge of the simulation box. 

All calculations start from a perfect fee crystal. Our runs consisted of a total of about 
5.5 • 10 5 passes during each of which an attempt was made to move particles and shadows. 
We allowed about 50 • 10 3 passes for equilibration followed by about 5 • 10 5 passes which 
comprise the equilibrated random walk. 

In Table | we show the energy per particle obtained from the M7+AS shadow wave- 
function after simulations with N = 108 particles at several densities of liquid 4 He. Also 
included in the table are several results from the literature, GFMC refers to the results of 
Kalos et affl. 

In Table [TT] we show the values of the optimum variational parameters b, C, r, and a for 
the M7+AS shadow wavefunction at different densities p in the liquid phase. 

The energy per particle for the M7+AS shadow wavefunction after simulations with 



N = 108 particles at some densities in the solid phase is shown in Table |IJ. In the same 
table we show the VMC results obtained with the M+MS and M+AS shadow wavefunction, 
as well as the GFMC results. 



Table [IV| shows the values of the optimum variational parameters b, C, r, and a for the 



M7+AS shadow wavefunction at different densities p in the solid phase. We fit our equation 



of state in the liquid phase to a cubic polynomial of the form 



E(p) = E + B 



P~ Po 
Po 
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C 



P~ Po 
Po 
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(11) 



where po is the equilibrium density. A similar function has been used to fit the experimental 
equation of stateE£l and we use it to analyze our M7+AS results. The values of the param- 
eters in the fit and their errors are shown in Table [V], together with other results from the 
literature. One notes that the values of the coefficients B, C, and p^a 3 for the M7+AS 
shadow wavefunction are in good agreement with both GFMC and experimental results!^. 
In the solid phase we used the same parametrization as reported for the M+AS case. We 
fited the energy to a cubic polynomial of the form 



E{p) =E + B 



P~ Ps 



Ps 



2 



+ c 



P- Ps 
Ps 



:-! 



(12) 



where the specific density pg<r 3 =0.4486 is taken from the GFMC calculation. The values of 



the parameters to this fit are shown in Table VI, together with other results from the liter- 



ature. Again the M7+AS shadow wavefunction shows a good agreement with GFMC with 
the exception of a discrepancy in the coefficient C. In Fig [J] we plot the equation of state for 
4 He liquid as obtained by using the values of the fitting parameters reported in Table |V|. As 
already known, one notes that both M+AS and M7+AS wavefunctions give a better equa- 
tion of state than the shadow wavefunction with fully optimized Jastrow particle-particle 
correlations (OJ+AS), although the OJ+AS wavefunction gives somewhat lower energies. 
It has been argued^ that possible causes of such behavior are the incomplete determination 
of the coefficients in the basis-set expansion for the OJ+AS wavefunction, or the missing 
full reoptimization of the shadow parameters. Indeed the recent VMC calculations with a 
fully optimized shadow wavefunctionliZI confirm this latter possibility. 

In Fig. Q, we show the radial distribution function g(r) obtained at the GFMC equilibrium 
density pa 3 = 0.365. Our maximum of g(r) is obtained at the same position r max as the 
GFMC value and it is clear tha there is no shifting of our curve to larger values of r. Our 
variational peak g(r max ) is a little smaller than the GFMC value. Fig. |3| shows our results for 
the radial distribution function g{r) determined at the GFMC freezing density pa 3 = 0.438. 
Statistical errors in the GFMC g(r) near the maximum g(r max ) are large, so a detailed 
comparison with GFMC is not possible. It appears that the position of our maximum of 



g(r) compares very well with the GFMC value, but the variational peak g(r max ) is again 
smaller. The trend seen at the GFMC equilibrium and freezing densities is repeated at all 
other densities. 

In Fig. |] we show S(k) at the equilibrium density pa 3 = 0.365. The experimental 
S(k) shown in this figure is the result reported by Svensson et al^ obtained by neutron 
diffraction at saturated vapor pressure at T = 1.0K°. The agreement of the variational 
structure factor with experiment is seen to be very good for all k-s, except for small k. This 
is to be expected since our M7+AS wavefunction does not contain the proper long-range 
correlations necessary for the linear behavioral of S(k) which is observed in 4 He. 

In this work we have demonstrated the utility of our M7+AS shadow wavefunction for 
studying the ground state properties of liquid and solid 4 He. The use of an inverse seventh 
power as the particle-particle correlation function has significantly improved the ground 
state energy, the equation of state and the radial density distribution. Such an improvement 
was obtained with very little additional computational effort. The wavefunction remained 
simple, compact and portable. 
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TABLES 
TABLE I. Energies in liquid 4 He at several densities including the experimental equilibrium 
density (pa 3 = 0.365) and the GFMC freezing density (pa 3 = 0.438) with a = 2.5561. All simula- 
tions use the Aziz HFDHE2 potential and have been performed for systems of N = 108 particles. 
The energies are given in Kelvin per particle. M+MS refers to a shadow wavefunctiorJHI with McMil- 
lan fifth power-law pseudopotential (m=5) for both particle-particle and shadow-shadow pseudopo- 
tentials. M+AS refers to a shadow wavefunctionE3 with a rescaled Aziz HFDHE2 shadow-shadow 
pseudopotential and a McMillan fifth power-law particle-particle pseudopotential (m=5). M7+AS 
refers to a shadow wavefunction with a rescaled Aziz HFDHE2 shadow-shadow pseudopotential 
and a McMillan seventh power-law particle-particle pseudopotential (m=7) as used in this work. 
GFMC refers to the Green's Function Monte Carlo calculationsH with the Mcmillan fifth power-law 
form for the importance and starting function. 



pa 3 


Method 


Trial function 


Energy (K°) 


0.328 


VMC 


M+AS 


-6.561 ± 0.032 




VMC 


M7+AS 


-6.571 ± 0.015 




GFMC 




-7.034 ± 0.037 


0.365 


VMC 


M+MS 


-6.061 ± 0.025 




VMC 


M+AS 


-6.599 ± 0.034 




VMC 


M7+AS 


-6.664 ± 0.021 




GFMC 




-7.120 ± 0.024 


0.401 


VMC 


M+AS 


-6.398 ± 0.019 




VMC 


M7+AS 


-6.497 ± 0.012 




GFMC 




-6.894 ± 0.048 


0.438 


VMC 


M+MS 


-5.360 ± 0.035 




VMC 


M+AS 


-5.871 ± 0.016 




VMC 


M7+AS 


-6.067 ± 0.010 




GFMC 




-6.564 ± 0.058 
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TABLE II. The variational parameters of the M7+AS shadow wavefunction used in the sim- 
ulation of 4 He liquid with N = 108 particles at different densities. 



pa" 


b/a 


Ca s 


r(K-i) 


a 


0.328 


1.02 


5.5 


0.088 


0.915 


0.365 


1.01 


5.5 


0.095 


0.915 


0.401 


1.01 


6.0 


0.105 


0.920 


0.438 


1.01 


5.9 


0.110 


0.910 
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TABLE III. Energies in solid 4 He at several densities including the GFMC melting density 
(pa 3 = 0.491) with a = 2.556A All simulations use the Aziz HFDHE2 potential and have been 
performed for systems of A = 108 particles. The energies are given in Kelvin per particle. The 
notation is the same as in Table |. The GFMC result at density pa 3 = 0.550 was interpolated from 
the GFMC results at pa 3 = 0.526 and pa 3 = 0.560. The M7+AS results represent this work. 



pa 3 


Method 


Trial function 


Energy (K°) 


0.491 


VMC 


M+MS 


-5.004 ± 0.055 




VMC 


M+AS 


-5.052 ± 0.014 




VMC 


M7+AS 


-5.324 ± 0.010 




GFMC 




-5.610 ± 0.030 


0.550 


VMC 


M+MS 


-3.521 ± 0.032 




VMC 


M+AS 


-3.639 ± 0.012 




VMC 


M7+AS 


-3.724 ± 0.017 




GFMC 




-4.197 ± 0.030 


0.589 


VMC 


M+AS 


-1.947 ± 0.012 




VMC 


M7+AS 


-2.097 ± 0.010 




GFMC 




-2.680 ± 0.060 



TABLE IV. The variational parameters of the M7+AS shadow wavefunction used in the 
simulation of 4 He solid with N = 108 particles at different densities. 



pa 3 


b/a 


Ca 3 


r(K^) 


a 


0.491 


1.00 


5.7 


0.110 


0.875 


0.550 


1.00 


5.9 


0.100 


0.890 


0.589 


1.00 


6.5 


0.110 


0.900 
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TABLE V. Fit parameters of the equation of state for 4 He in the liquid phase. The OJ+AS 
shadow wavefunctionEj incorporates a fully optimized Jastrow particle-particle pseudopotential. 
The GFMC result is taken from Kalos at alQ. The experimental equation of state (Exp) is taken 
from Roach et ali3. 





E (K°) 


B(K°) 


C(K°) 


Po^ 3 


M+AS 


-6.610 ± 0.036 


10.3 ± 5.5 


11.3 ± 18.5 


0.3535 ± 0.0043 


OJ+AS 


-6.796 ± 0.025 


14.10 ± 4.18 


-18.7 ± 18.1 


0.3567 ± 0.0032 


M7+AS 


-6.662 ± 0.020 


14.08 ± 1.36 


6.72 ± 8.02 


0.361 ± 0.001 


GFMC 


-7.110 ± 0.023 


10.08 ± 3.20 


12.59 ± 8.50 


0.3600 ± 0.0049 


Exp 


-7.14 


13.65 


7.67 


0.365 



TABLE VI. Fit parameters of the equation of state for 4 He in the solid phase. The notation 
is the same as in Table [V|. The fitting curve is E = Eq + B[(p — p s )/p s ] 2 + C[(p — p s )/p s ] 3 , where 
p s a 3 = 0.4486 is taken from the GFMC result. 





E (K°) 


B(K°) 


C(K°) 


Ps° 3 


M+AS 


-5.340 ± 0.021 


31.00 ± 1.50 


9.92 ± 4.34 


0.4486 ± 0.0097 


OJ+AS 


-5.81 ± 0.02 


47.7 ± 1.4 


-33.89 ± 4.14 


0.4486 ± 0.0097 


M7+AS 


-5.693 ± 0.013 


32.92 ± 1.00 


-13.92 ± 1.45 


0.4486 ± 0.0097 


GFMC 


-5.899 ± 0.121 


31.95 ± 5.26 


3.395 ± 80.0 


0.4486 ± 0.0097 
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FIGURES 
FIG. 1. Equation of state of 4 He liquid. The lines are polynomial cubic fits. The solid line is 
the experimental resultslij. The dotted line denotes the GFMC results. The dashed line denotes 
OJ+AScJ. The solid line with opaque circles is the M7+AS result. The solid line with crosses is 
the M+AS resultEl. 

FIG. 2. Radial distribution function g(r) at the GFMC equilibrium density pa 3 = 0.365 after 
a VMC simulation with N=108 particles. The solid line denotes Ml + AS. Filled circles are the 
GFMC results of Kalos et aM. 

FIG. 3. Radial distribution function g{r) at the GFMC freezing density pa 3 = 0.438 after a 
VMC simulation with N=108 particles. The solid line denotes M7 + AS. Filled circles are the 
GFMC results of Kalos et aM. 

FIG. 4. Static structure factor S(k) of liquid 4 He at the GFMC equilibrium density pa 3 = 0.365. 
The filled circles show our results obtained from the formula S(k) = ^(/O-kPk) at a discrete set of 
/c-points. The solid line denotes the experimental results reported by Svensson and co-workeral3 
obtained at saturated vapor pressure by means of neutron diffraction at temperature T = 1.0K°. 
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